#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <time.h>
#define ITERACOES 25
#define DEZMI 10000000

typedef struct{
    mpf_t a;
    mpf_t b;
    mpf_t t;
    mpf_t p;
}Data;

int main(){
    time_t start, end;
    Data D[ITERACOES];
    int i;
    int n;

    mpf_set_default_prec((DEZMI*4)+16);
    mpf_t PI;
    mpf_t temp;

    //start = time(NULL);

    for(i=0; i<ITERACOES; i++){
        //mpf_init(D[i].a);
        mpf_init_set_d(D[i].a, 0.0);

        //mpf_init(D[i].b);
        mpf_init_set_d(D[i].b, 0.0);

        //mpf_init(D[i].t);
        mpf_init_set_d(D[i].t, 0.0);

        //mpf_init(D[i].p);
        mpf_init_set_d(D[i].p, 0.0);
    }

    //mpf_init(PI);
    //mpf_init(temp);
    mpf_init_set_ui(PI, 0L); // pi =0
    mpf_init_set_ui(temp, 0L); //temp = 0


    mpf_init_set_ui(D[0].a, 1L);

    mpf_sqrt_ui(temp, 2L); // temp = sqrt(2)
    mpf_ui_div(D[0].b, 1L, temp); // b0 = 1/sqrt(2)

    mpf_init_set_d(D[0].t, 0.25); //t0 = 0.25

    mpf_init_set_ui(D[0].p, 1L); // p0 = 1

    /*And let's the magic begin: */

    for(n=0; n<ITERACOES-1; n++){
        //Value for A(n+1)
        mpf_add(temp, D[n].a, D[n].b);
        mpf_div_ui(D[n+1].a, temp, 2L);

        //Value for B(n+1)
        mpf_mul(temp, D[n].a, D[n].b);
        mpf_sqrt(D[n+1].b, temp);

        //value for T(n+1)
        mpf_sub(temp, D[n].a, D[n+1].a);
        mpf_pow_ui(temp, temp, 2L);
        mpf_mul(temp, D[n].p, temp);
        mpf_sub(D[n+1].t, D[n].t, temp);

        //value for P(n+1)
        mpf_mul_ui(D[n+1].p, D[n].p, 2L);

        mpf_clear(D[n].a);
        mpf_clear(D[n].b);
        mpf_clear(D[n].p);
        mpf_clear(D[n].t);
    }

    //value for PI
    mpf_add(temp, D[n].a, D[n].b);
    mpf_pow_ui(PI, temp, 2L);
    mpf_mul_ui(temp, D[n].t, 4L);
    mpf_div(PI, PI, temp);

    //end = time(NULL);
    //printf("Execution Time: %.10fs\n", difftime(end,start));

    FILE *arquivo;
    arquivo=fopen ("out_Gauss.txt","w");
    if (!arquivo)
        printf ("Erro na abertura do arquivo de saida com valor de pi.");
    
	mpf_out_str( arquivo, 10, (DEZMI*4)+16, PI );
    fclose(arquivo);

    //gmp_printf("PI Value Stored in file, about : %.10Ff \n", PI);


    return 0;
}
